function [tOut, uOut, udtOut, udtdtOut] = h2hydro(params, uIn)



  % ****************************
  % run the ode solver
  % ****************************

  time_steps = (0:params.Tsteps-1)/params.Tsteps*params.Tfinal;


  tic;
  %uOut   = zeros(params.Tsteps, Nt.*5);
  %tOut   = zeros(params.Tsteps, Nt.*5);

    opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'OutputFcn', @(t,u,flag) hplot(t, u, flag, params));


  switch params.hydroOrder
    case 2
      eom = @heom2;
    case -2
      eom = @heom2nr;
    case 0
      eom = @heom0;
  end

  [tOut, uOut] = ode45(@(t,u) eom(t,u,params), time_steps, uIn, opts);

  elapsed = toc;
  disp(sprintf('computed hydro simulation in %fsec', elapsed));


  % ****************************
  % compute 1st time derivatives
  % ****************************

  tic;

  Nt = params.Nx*params.Ny;
  udtOut   = zeros(params.Tsteps, Nt.*5);

  % first derivative from eom
  for it=1:params.Tsteps
    t = tOut(it);
    u = uOut(it, :);
    
    udtOut(it,:) = eom(t, u, params);
  end

  % ****************************
  %  get rid of Pi and B
  % ****************************

  uOut  (:,3*Nt+1 : 5*Nt) = [];
  udtOut(:,3*Nt+1 : 5*Nt) = [];


  % ****************************
  % compute 2nd time derivatives
  % ****************************

  derorder = params.timeOrder;
  % second der from finite differences, O(8)
  derker = d1ker(params.timeOrder);
  %derker = [1.0/12; -2.0/3; 0; 2.0/3; -1.0/12]; % o4
  udtdtOut = conv2(udtOut, derker, 'same')./params.Dt;

  % delete derivative padding
  deli = [1:(derorder/2) (length(tOut)-derorder/2+1):length(tOut)];
  tOut    (deli, :) = [];
  uOut    (deli, :) = [];
  udtOut  (deli, :) = [];
  udtdtOut(deli, :) = [];
  elapsed = toc;
  disp(sprintf('computed higher hydro time derivatives in %fsec', elapsed));

  params.Tbegin = params.Tbegin + derorder./2.*params.Dt;
  params.Tfinal = params.Tfinal - derorder./2.*params.Dt;
  params.Tsteps = params.Tsteps - derorder;


  % ****************************
  % compute ut
  % ****************************

  ux = uOut(:, 1*Nt+1 : 2*Nt);
  uy = uOut(:, 2*Nt+1 : 3*Nt);

  uxdt = udtOut(:, 1*Nt+1 : 2*Nt);
  uydt = udtOut(:, 2*Nt+1 : 3*Nt);

  uxdtdt = udtdtOut(:, 1*Nt+1 : 2*Nt);
  uydtdt = udtdtOut(:, 2*Nt+1 : 3*Nt);

  ut = sqrt(1+ux.^2+uy.^2);
  utdt = (uxdt.*ux + uydt.*uy)./ut;
  utdtdt = (1+ux.^2+uy.^2).^(-3/2).*(uxdt.^2.*(1+uy.^2)+(-2).*ux.*uxdt.*uy.*uydt+(1+ux.^2).*uydt.^2+(1+ux.^2+uy.^2).*(ux.*uxdtdt+uy.*uydtdt));

  uOut    (:, 3*Nt+1 : 4*Nt) = ut;
  udtOut  (:, 3*Nt+1 : 4*Nt) = utdt;
  udtdtOut(:, 3*Nt+1 : 4*Nt) = utdtdt;



end
